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Abstract. Computing spherical liarmonic decompositions is a ubiquitous tectmique that arises 
in a wide variety of disciplines and a large number of scientific codes. Because spherical harmonics are 
defined by integrals over spheres, however, one must perform some sort of interpolation in order to 
compute them when data is stored on a cubic lattice. Misner (2004, Class. Quant. Grav., 21, S243) 
presented a novel algorithm for computing the spherical harmonic components of data represented on 
a cubic grid, which has been found in real applications to be both efficient and robust to the presence 
of mesh refinement boundaries. At the same time, however, practical applications of the algorithm 
require knowledge of how the truncation errors of the algorithm depend on the various parameters 
in the algorithm. Based on analytic arguments and experience using the algorithm in real numerical 
simulations, I explore these dependencies and provide a rule of thumb for choosing the parameters 
based on the truncation errors of the underlying data. I also demonstrate that symmetries in the 
spherical harmonics themselves allow for an even more efficient implementation of the algorithm than 
was suggested by Misner in his original paper. 

1. Introduction. Spherical harmonic decomposition is a ubiquitous mathemat- 
ical technique that arises in a wide variety of disciplines. It is no surprise, therefore, 
that it also arises in an equally varied range of scientific computing applications. It 
plays a key role, for example, in currently active algorithms, for determining com- 
patible docking configurations of organic molecules [S] 112) , modeling atmospheric and 
oceanic flows ^ |51 lll| , representing human population density on the surface of the 
earth |1()| . analyzing experimental data regarding turbulent magnetohydrodynamic 
flows "H^ , and analyzing numerical simulations of gravitational radiation emitted from 
the collisions of black holes |21IS1- 

In numerical calculations structured on a cubic grid, however, extracting spherical 
harmonic components can be non-trivial since the spherical harmonic components 
of a function $ are defined by integrals over spheres 



that have few if any intersections with the cubic lattice. Misner recently presented 
a particularly nice approach to solving this problem, which does not require explicit 
interpolations from the cubic grid to an integration sphere The mathematical 
grounding of the algorithm is laid out completely in Ref. [S], but the original refer- 
ence does not provide any detailed analysis of the numerical errors incurred by the 
approximation. In practical applications, such error estimates have been found essen- 
tial for choosing the parameters of the algorithm and for understanding the results 
of simulations In addition, making use of the symmetries of the spherical har- 

monics allows certain simplifications of the algorithm as applied to generic data and, 
especially, as applied to data with explicit grid symmetries when only the independent 
portion of the data is evolved. 

It should be noted that although, in my discussion, I will always speak of spherical 
harmonics in their usual sense, everything in this paper is also true for spin-weighted 
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spherical harmonics (see Refs. 0I[7| for definitions) except for Section[Sl spin- weighted 
spherical harmonics do not generally have the simple symmetries that the (usual) 
scalar spherical harmonics have. Note also that I follow Misner in treating only 
the case of uniform grids, although this method was successfully applied to non- 
uniform grids using fixed mesh refinement. For more details on applications using 
spin- weighted harmonics and run on non- uniform grids, see Refs. !3|II3]- 

The rest of the paper is organized as follows: In Section |21 I provide a brief 
summary of Misner algorithm in order to establish notation for the following sections. 
The error analysis for the algorithm is contained in Section |21 which leads directly 
to Sectional where I lay out a rule of thumb for choosing the various parameters in 
the algorithm. In Section |S1 I provide a detailed analysis of the issues associated with 
and simplifications resulting from various symmetries in the spherical harmonics, and 
how those symmetries can be exploited for data that possesses explicit symmetries. 

2. Methodology. In this section I outline in some detail the steps of the Misner 
algorithm, although I do not provide a full derivation. (See Ref. [HI for the derivation.) 
The purpose of this chapter is primarily to establish notation and conventions for the 
rest of the paper. 

In order to begin, two definitions are need. First define a radial function 

i?„(r;i?,A) =r-iy^P„ (2.1) 

in terms of the usual Legendre polynomials P„. Here R and A are parameters that 
will be associated with the radius at which the spherical harmonic decomposition is 
desired and half of the thickness of a shell centered on that radius. From this, define 

Ynl„^{r,e,(j)) = Rn{r)Yi^{e,(j)) (2.2) 

which form a complete, orthonormal set with respect to the inner product 

(/l.9> = / K^)9{x)d?x (2.3) 

JS 

on the shell S — {(?', 0,4)) \ r e [R — A, R + A]}. Note also that, because the functions 
Rn form a complete set on the shell, 

YU0,cl>)<!>{r,9,4)d''x (2.4) 
and that the term in brackets 

oo 

Rn{R\ R, A)Rn{r; R, A) = r-^5{R - r) (2.5) 

n=0 

is a delta function.^ (Compare (|2.4|l to (|l.l|l .') 

On a finite grid F, the inner product (|2.3|l will have the form 

^The Dirac delta function is defined by the property that, for any function /, the integral 
Ja fiy)^{^ ~ y)d-y f{x) when x G (a, b) and is otherwise. 
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^lm{t,R)= J 



n=0 



Rn{R;R, A)Rn{r;R, A) 




where each point has some weight w^- This weight was given the form 

\r-R\> A + h/2 
\r~R\<A^h/2 (2.7) 
\R ~ r\)h'^ otherwise 

by Misner, where h is the grid spacing. Only cases with A > h/2 are considered. This 
means, roughly, that points entirely within the shell S are weighted by their finite 
volume on the numerical grid, points entirely outside of the shell S have zero weight, 
and points near the boundary are weighted according to the fraction of their volume 
inside S."^ 

With the numerical inner product H2.6|l . and letting capital Roman letters A = 
(nlm) represent index groups, the Ya are no longer orthonormal. Their inner product 



{Ya\Yb) = Gab = Gba (2.8) 

forms a metric for functions on the shell. Although a priori this matrix appears to 
be complex valued, I will show in Section |S1 that it is actually real-symmetric and 
sparse. For now it suffices to follow Misner in denoting it as generically Hermitian. 
The inverse to this metric G^^ can be used to raise indices on functions defined on 
the sphere. 

Making use of this new metric, and with some further analysis, the approximation 
for the spherical harmonic coefficients 

$,™(t, R)^Y. Rim{x] R)wx^{t, x) (2.9) 

follows with 

N 

i?z,„(x;i?) = ^i?„(i?)r"'™(a;) (2.10) 

in terms of Y^ = G^'^Yb, not Ya. 

3. Error Analysis. In Ref. [S], Misner provides an algorithm for computing 
spherical harmonic components on a cubic lattice, but he does not provide any analysis 
of the errors involved. In a practical application of the method, such error analysis, 
in particular an analysis of the convergence order in grid spacing, is very valuable. 
For this reason, I wish to examine the issue more closely here. 

To begin the analysis, note that there are three parameters in the algorithm, the 
grid spacing h, the width of the shell A, and the number of terms in the sum over 
basis polynomials N. (Compare H2.4|) to (|2.9|) and H2.10|) .'l In the limit that /i — > and 
iV — > oo, the numerical algorithm goes to the continuum theory. Note that, when h 
and N go to their limits, any value for A is allowed, since l|2.4|l is an exact, continuum 
expression. For a finite value of N, however, the quality of the approximation in the 
radial direction is a function of both N and A. Moreover, tying the value of A to the 
grid spacing /i, though not necessary from fundamental considerations, does in fact 
have advantages in terms of convergence behavior that will become clear below. 

^Actually the boundary points are weighted by the fraction of their volume that would be inside 
S if the point were on a coordinate axis. See Ref. for more details on the definition of the weights. 
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Table 3.1 



The first few non-trivial values of cjv fc(l), which are the coefficients of the Taylor expansion 
of a function integrated against d{x; N, 0, 1). (In general, cpf fc(A) = cjv fc(l)A'° .) The values of k 
run across and the values of N run down. The fact that the first coefficient is always 1, and that, 
by increasing N, more of the sub-leading coefficients are indicates that increasing N increases the 
order of convergence of the Misner algorithm ( provided that A oc h). 



I first consider the behavior of the algorithm as a function of A for fixed values 
of A^. Define 

v-^ 2n+l f X — R\ , , , , 

d(x;^,i?,A) = ^^^Pj^- P„(0) (3.1) 

n=0 ^ 

which is closely related to the delta function H2.5() in the limit N ^ oo. For finite N, 
this should approximate the delta function to some order of accuracy. To quantify 
this, consider the approximation 

/(O) « InU] = r d{x; N, 0, A)fix)dx (3.2) 

for a suitably smooth test function /. In the limit that N oo, this is exact. For 
any finite N, this expression can be analyzed by Taylor expanding the test function 
around the point a; = 0. This gives, after rearranging terms and noting that all of the 
odd terms vanish by symmetry, 

oo 

k=0 ^ 

where 

.A 

CN,kiA)^ x''d{x;N,0,A)dx (3.4) 

J-A 

are the coefficients of the expansion. In order to have a high order method, I need 
Cn,o = 1; a-iid the coefficients corresponding to the next few values of k to vanish. 
Table im shows the first few coefficients CN,k{^) = A^'^cjv.fc(A) for even N in [0,10]. It 
is clear that each time N is increased, one more of the sub-leading coefficients vanishes. 
This means that there are error terms proportional to A, and that the leading order 
term arises at 0(A^+^). In order to prevent this term, at high resolution, from 
dominating over discretization errors scaling like powers of the grid spacing, one should 
choose A cx ft,. 

Additional error terms proportional the grid spacing h depend primarily on the 
order of accuracy in the volume integral when using the weights defined by (|2.7|) . For 

4 



a finite volume, this weighing scheme provides an approximation that scales as 0{h), 
but, for a region that is also scaling with h, the resulting integral scales as 0{h?). 
This provides additional motivation for choosing A oc /i since most applications will 
require at least second order accuracy. 

For higher than second order accuracy, a new scheme for computing the volume 
integrals is needed, but the rest of the analysis here holds true. Given such a scheme, 
the analysis here shows how to choose the remaining parameters to ensure that nu- 
merical errors scale like any desired power of the grid spacing. 

4. Choosing the parameters. In practice, the grid spacing parameter h is 
usually chosen to resolve the sources without exceeding the physical limits of the 
computer. I would not expect, in general, that the grid spacing would be chosen 
based on the needs of this algorithm. For that reason, let me assume now that h is 
chosen, and discuss how to choose the remaining parameters N and A. In this section 
I will discuss some of the theoretical issues that should be considered when choosing 
the parameters, leading to a rule of thumb that is valid based on this analysis and my 
experience with the algorithm. 

The error analysis of Section 13 implies that for fixed A, increasing the value of 
N decreases the error term. It also implies that for fixed N , increasing the value of 
A increases the error term. This suggests taking A as small as possible, and N as 
large as possible to make the error term as small as possible. This must be balanced, 
however, against practical limitations. Certainly the shell thickness A needs to be 
large enough so that there are some grid points within the shell, otherwise the whole 
procedure is undefined. For fixed iV, a stronger restriction requires that the Legendre 
polynomial Pm can be resolved over the shell. Without this condition, there would 
seem to be no benefit to taking higher values of N . Getting higher accuracy in practice 
requires finding a proper balance between choosing A small and N large. 

In making this balance, however, one must keep in mind that the error in the 
method is partially determined by the weighting scheme (|2.7|l . which is only second 
order accurate in the grid spacing. I am, in addition, going to choose A cx ft, for reasons 
described above. This already suggests that taking N larger than two is pointless, 
since choosing N — 2 already makes the piece of the error that is proportional to A 
scale like 0{A*) (cf. Table IXT|) . meaning that it will be an error term of sub-leading 
order in grid spacing. But once this term is of sub-leading order, it is much less 
important how large I choose A, provided that I still choose it proportional to the 
grid spacing. I therefore adopt the following 

Rule of Thumb: Choose N just large enough to ensure that the 
error term proportional to A is an error term of sub-leading order in 
grid spacing. Choose A just large enough to safely resolve Pn on the 
shell. 

With this rule of thumb, and the second order accurate weighting scheme (|2.7|l . I 
found the choices N — 2 and A = 3h/A completely satisfactory for a second order 
accurate code. Note that this corresponds to Misner's choice of A in Ref. With 
= 2 I found that larger values of A are also acceptable. Numerical results justifying 
these estimates appear in Ref. 

This is illustrated in Figure 14.11 which shows the errors in a test function as 
compared to an analytic solution. (The physical problem is the propagation of a 
linear gravitational wave through a mesh refined grid, as described in Refs. It 
is clear from the figure that choosing N large enough makes a dramatic impact on 
the truncation error, whereas the exact value of A for fixed N is less important. In 
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Comparision of Error in Extracted Wave Using Misner's Ivlethod witli Various Parameters (r=3 X) 
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Fig. 4.1. A comparison of numerical results as a function of parameters. The lines show errors 
in a numerical solution as compared to an analytic solution for a sample problem. When N = 2 the 
errors are smaller than when N = 0, consistent with the fact that the errors scale like the square of 
grid spacing for N = 2 and only like the grid spacing for N = 0. For fixed N, the size of A is fairly 
unimportant. 



this test case, the underlying simulation was only second order accurate, so there is 
no benefit from choosing N > 2. 



5. Symmetry issues. There are two points of interest related to this method of 
spherical harmonic decomposition and symmetries. The first was mentioned briefly in 
Section 12 namely that symmetry causes the metric Gab to be real and sparse. The 
second deals with implementing the method for grids in which explicit symmetries are 
enforced on grid functions in order to reduce the computational load of the simulation. 
In these cases, in which data is not evolved over a whole extraction sphere, additional 
analysis is required to demonstrate that the method is well defined and to understand 
how to most efficiently implement it. The primary result on this second topic is that 
the adjoint harmonics of (|2.1()(l have the same symmetries under reflection as the 
original spherical harmonics Ya- 

Consider first the implications of symmetry on the metric Gab- The symmetries 
of the spherical harmonics, summarized in Table 15.11 cause the imaginary part of 
all terms in the integral 1)2.8(1 to cancel in pairs of points on the sphere related by 
reflections through coordinate planes. The reason is that each of the four signs (+1, 
(— 1)™, (—1)', and (—1)'+™) appears twice in Tablc lSlTl once for a term that is complex 
conjugated and once for a term that is not. The matrix is also sparse. By similar 
reasoning, for certain values of I and m, the terms in the integral ((2.8(1 can cancel in 
sets of four. Both of these facts can be seen at once through a simple calculation. 
The idea is to break the integral into parts using the second and third columns of 
Table I^TI and then to simplify using the last two columns. Considering first just the 
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Table 5.1 



The table shows how the arguments of spherical harmonics transform under reflections through 
various Cartesian planes. The first column indicates which coordinates have their signs inverted, 
while the second and third columns give the new angular arguments to the spherical harmonic Yi^j^. 
Alternatively, the fourth and fifth column show, respectively, the overall sign in front of and whether 
or not to complex conjugate the given spherical harmonic with the original angular arguments. The 
second row, for example, says that Yi^{—x,y,z) = Yi^{9,it — <f>) = (6, 0), where (9,0) 

are the angular coordinates of the point {x,y,z). Note that this table differs slightly from that in 
Ref. I^. The table here is correct. 



symmetries under reflection through the xy-plane and recaUing the definition H2.8|) , 

GaB^ j>%mAS.4>)yi2rnAG,^)dn (5.1a) 

Yi,mA&,4,)Yi,^,{9,c^)dn 

/ Yi„nA7^-9,C^)Yi,„,,{TT-9,C^)dn (5.1b) 

Jo 

n2w nw/2 

= [l+{-iy^+'-] / Yl,„^A0,m2mA^A)dn. (5.1c) 

Jo Jo 

(I have suppressed the radial functions since they play no role here.) Repeating the 
procedure for reflections through the xz- and yz-planes shows that 

Gab — 2(Tmi+m2,ii+i2 

/ / Re{Yi,.rm{9,m2mA0,^)}dn (5.2) 
Jo Jo 

where 

'^mMM+h = 1 + (-1)™^+™^ + (-1)'^+'^ + (_l)™i+™2+h+i2. (5.3) 

This proves that the matrix is real. In addition, the matrix element is zero by sym- 
metry whenever 

0'mi+m2,il+i2 = (5-4) 

which is true for 56 of the 81 matrix elements that exist when considering a fixed 
value of n and all values of / and m for I < 2. Of the remaining 25 matrix elements, 
9, of course, are the diagonal elements that go to unity in the continuum limit. The 
exact break-down of which such elements must be zero by symmetry is summarized 
in Table Knowing that the matrix is real-symmetric and sparse allows for a more 
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Table 5.2 



The table summarizes which entries of Gab identically vanish because of the symmetries of the 
spherical harmonics under reflections through coordinate planes for all values of I and m with I < 2. 
(This is governed by equation (5.4)-) Of the 81 possible matrix elements, only 25 have non-trivial 
values. 



efficient implementation of the algoritlim in general. It is also extremely useful in 
analyzing the algorithm in the context of the second topic of this section, explicit grid 
symmetries. 

When evolving initial data with known symmetries, it is very common to evolve 
only that part of the data that is unique. In such cases, an appropriate symmetry 
boundary condition is applied at some edges of the grid. This is, however, inconvenient 
for wave extraction since computing spherical harmonic components (by any method) 
requires integrating over the full sphere. If data with octant symmetry, for example, 
is evolved only in a single octant, it is neither sufficient to apply the decomposition 
algorithm in that one octant nor to multiply the result of a single octant by 8 since 
the symmetry may forbid some modes as well as repeating them. 

In principle the problem appears to be even more difficult for this particular de- 
composition method. Although the spherical harmonics have well defined symmetries 
under reflections, as summarized in Table [FTl it is the adjoint harmonics that appear 
in (j2.10|l . The adjoint harmonics, however, are constructed by contracting G^^ with 
the usual spherical harmonics, and this appears to mix different values of I and m. 
While this mixing does occur, the the matrix Gab is sparse in just the right way to 
ensure that the adjoint harmonics have the same symmetries at the usual spherical 
harmonics. 

A particular choice of the mapping (n, Z, m) h-s- A makes this easiest to see. Specif- 
ically, considering all values of / and m with I <2, there is a basis in which Gab takes 
block diagonal form 



(Gab) 



( Si 



V 



S4 } 



(5.5) 



with all unwritten entries identically zero by symmetry. In this expression Si is a 
4iV X 4iV matrix over the basis functions B\ — {l^noo, ^n2.-2, i^ri20i i^n22}; S2 is an 
N X N matrix over the basis functions B2 = {i^nio}; S3 is a 2A^ x 2N matrix over 
the basis functions B3 — {y„i^_i, K^n}; and S4 is a 2N x 2A^ matrix over the basis 
functions B4 = {i^n2,-i, i^n2i}- In this basis the matrix is block diagonal, so the 
inverse matrix 



1 



(G 



AB\ 



— -1 
"3 



sr^ / 



(5.6) 



is also block diagonal and the different basis sectors do not mix. This last point is key. 
It implies that any particular adjoint harmonic is a linear combination of spherical 
harmonics from a single set Bk 



where k is the index such that Ynim G Bk- Because, in each set Bk, the parity 
of I and the parity of m is the same on each Ynim G Bk, and because, in light of 
Table 15.11 it is the parity of I and m that determines the symmetries of Ynim under 
reflections through planes, every spherical harmonic in Bk for any fixed k has the same 
symmetries under reflections as any other spherical harmonic in Bk- This implies that 
the adjoint harmonics also share this symmetry under reflection. 

6. Discussion. In this paper I have provided detailed error estimates of the 
Misner algorithm for computing spherical harmonic components of data represented 
on a cubic grid. This analysis allows one, in principle, to chose the parameters of the 
algorithm such that its numerical errors scale any desired power of the grid spacing. 
The only limitation of this in practice is finding a scheme for approximating volume 
integrals on a shell of sufficiently high accuracy. (Misner's original choice allows for a 
second order accurate result.) 

In addition, analysis of the symmetry properties of the spherical harmonics pro- 
vides two insights: First, the number of operations required to initialize the data 
structures used to compute spherical harmonic components can be reduced by com- 
puting only those elements of Gab that are not forbidden by symmetries, and, second, 
that the adjoint harmonics used by the algorithm have the same symmetries under 
reflections as the usual spherical harmonics. This second fact allows the method to 
be efficiently used on data with explicit grid symmetries when only the independent 
portion of that data is evolved. 
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